{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Modeling and Simulation in Python\n",
    "\n",
    "Chapter 22\n",
    "\n",
    "Copyright 2017 Allen Downey\n",
    "\n",
    "License: [Creative Commons Attribution 4.0 International](https://creativecommons.org/licenses/by/4.0)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Configure Jupyter so figures appear in the notebook\n",
    "%matplotlib inline\n",
    "\n",
    "# Configure Jupyter to display the assigned value after an assignment\n",
    "%config InteractiveShell.ast_node_interactivity='last_expr_or_assign'\n",
    "\n",
    "# import functions from the modsim.py module\n",
    "from modsim import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Vectors"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `Vector` object represents a vector quantity.  In the context of mechanics, vector quantities include position, velocity, acceleration, and force, all of which might be in 2D or 3D.\n",
    "\n",
    "You can define a `Vector` object without units, but if it represents a physical quantity, you will often want to attach units to it.\n",
    "\n",
    "I'll start by grabbing the units we'll need."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = UNITS.meter\n",
    "s = UNITS.second\n",
    "kg = UNITS.kilogram"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a two dimensional `Vector` in meters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "A = Vector(3, 4) * m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can access the elements by name."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.x"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.y"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The magnitude is the length of the vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.mag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The angle is the number of radians between the vector and the positive x axis."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we make another `Vector` with the same units,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "B = Vector(1, 2) * m"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can add `Vector` objects like this"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "A + B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And subtract like this:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "A - B"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can compute the Euclidean distance between two Vectors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.dist(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And the difference in angle"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.diff_angle(B)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we are given the magnitude and angle of a vector, what we have is the representation of the vector in polar coordinates."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "mag = A.mag\n",
    "angle = A.angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use `pol2cart` to convert from polar to Cartesian coordinates, and then use the Cartesian coordinates to make a `Vector` object.\n",
    "\n",
    "In this example, the `Vector` we get should have the same components as `A`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "x, y = pol2cart(angle, mag)\n",
    "Vector(x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Another way to represent the direction of `A` is a unit vector, which is a vector with magnitude 1 that points in the same direction as `A`.  You can compute a unit vector by dividing a vector by its magnitude:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "A / A.mag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Or by using the `hat` function, so named because unit vectors are conventionally decorated with a hat, like this: $\\hat{A}$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "A.hat()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Create a `Vector` named `a_grav` that represents acceleration due to gravity, with x component 0 and y component $-9.8$ meters / second$^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Degrees and radians\n",
    "\n",
    "Pint provides units to represent degree and radians."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "degree = UNITS.degree\n",
    "radian = UNITS.radian"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If you have an angle in degrees,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [],
   "source": [
    "angle = 45 * degree\n",
    "angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can convert to radians."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "angle_rad = angle.to(radian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If it's already in radians, `to` does the right thing."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "angle_rad.to(radian)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can also convert from radians to degrees."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "angle_rad.to(degree)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As an alterative, you can use `np.deg2rad`, which works with Pint quantities, but it also works with simple numbers and NumPy arrays:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "np.deg2rad(angle)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Create a `Vector` named `a_force` that represents acceleration due to a force of 0.5 Newton applied to an object with mass 0.3 kilograms, in a direction 45 degrees up from the positive x-axis.\n",
    "\n",
    "Add `a_force` to `a_grav` from the previous exercise.  If that addition succeeds, that means that the units are compatible.  Confirm that the total acceleration seems to make sense."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Baseball"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a `Params` object that contains parameters for the flight of a baseball."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "t_end = 10 * s\n",
    "dt = t_end / 100\n",
    "\n",
    "params = Params(x = 0 * m, \n",
    "                y = 1 * m,\n",
    "                g = 9.8 * m/s**2,\n",
    "                mass = 145e-3 * kg,\n",
    "                diameter = 73e-3 * m,\n",
    "                rho = 1.2 * kg/m**3,\n",
    "                C_d = 0.33,\n",
    "                angle = 45 * degree,\n",
    "                velocity = 40 * m / s,\n",
    "                t_end=t_end, dt=dt)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And here's the function that uses the `Params` object to make a `System` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [],
   "source": [
    "def make_system(params):\n",
    "    \"\"\"Make a system object.\n",
    "    \n",
    "    params: Params object with angle, velocity, x, y,\n",
    "               diameter, duration, g, mass, rho, and C_d\n",
    "               \n",
    "    returns: System object\n",
    "    \"\"\"\n",
    "    angle, velocity = params.angle, params.velocity\n",
    "    \n",
    "    # convert angle to degrees\n",
    "    theta = np.deg2rad(angle)\n",
    "    \n",
    "    # compute x and y components of velocity\n",
    "    vx, vy = pol2cart(theta, velocity)\n",
    "    \n",
    "    # make the initial state\n",
    "    R = Vector(params.x, params.y)\n",
    "    V = Vector(vx, vy)\n",
    "    init = State(R=R, V=V)\n",
    "    \n",
    "    # compute area from diameter\n",
    "    diameter = params.diameter\n",
    "    area = np.pi * (diameter/2)**2\n",
    "    \n",
    "    return System(params, init=init, area=area)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's how we use it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "system = make_system(params)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's a function that computes drag force using vectors:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "def drag_force(V, system):\n",
    "    \"\"\"Computes drag force in the opposite direction of `v`.\n",
    "    \n",
    "    V: velocity Vector\n",
    "    system: System object with rho, C_d, area\n",
    "    \n",
    "    returns: Vector drag force\n",
    "    \"\"\"\n",
    "    rho, C_d, area = system.rho, system.C_d, system.area\n",
    "    \n",
    "    mag = rho * V.mag**2 * C_d * area / 2\n",
    "    direction = -V.hat()\n",
    "    f_drag = direction * mag\n",
    "    return f_drag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can test it like this."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "V_test = Vector(10, 10) * m/s\n",
    "drag_force(V_test, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Here's the slope function that computes acceleration due to gravity and drag."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "def slope_func(state, t, system):\n",
    "    \"\"\"Computes derivatives of the state variables.\n",
    "    \n",
    "    state: State (x, y, x velocity, y velocity)\n",
    "    t: time\n",
    "    system: System object with g, rho, C_d, area, mass\n",
    "    \n",
    "    returns: sequence (vx, vy, ax, ay)\n",
    "    \"\"\"\n",
    "    R, V = state\n",
    "    mass, g = system.mass, system.g\n",
    "    \n",
    "    a_drag = drag_force(V, system) / mass\n",
    "    a_grav = Vector(0, -g)\n",
    "    \n",
    "    A = a_grav + a_drag\n",
    "    \n",
    "    return V, A"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Always test the slope function with the initial conditions."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "slope_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can use an event function to stop the simulation when the ball hits the ground:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "def event_func(state, t, system):\n",
    "    \"\"\"Stop when the y coordinate is 0.\n",
    "    \n",
    "    state: State object\n",
    "    t: time\n",
    "    system: System object\n",
    "    \n",
    "    returns: y coordinate\n",
    "    \"\"\"\n",
    "    R, V = state\n",
    "    return R.y"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "event_func(system.init, 0, system)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we can call `run_ode_solver`"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "results, details = run_ode_solver(system, slope_func, events=event_func)\n",
    "details"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final label tells us the flight time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [],
   "source": [
    "flight_time = get_last_label(results) * s"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The final value of `x` tells us the how far the ball landed from home plate:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [],
   "source": [
    "R_final = get_last_value(results.R)\n",
    "x_dist = R_final.x"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Visualizing the results\n",
    "\n",
    "The simplest way to visualize the results is to plot x and y as functions of time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 38,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = results.R.extract('x')\n",
    "ys = results.R.extract('y')\n",
    "\n",
    "xs.plot()\n",
    "ys.plot()\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Position (m)')\n",
    "\n",
    "savefig('figs/chap22-fig01.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We can plot the velocities the same way."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 39,
   "metadata": {},
   "outputs": [],
   "source": [
    "vx = results.V.extract('x')\n",
    "vy = results.V.extract('y')\n",
    "\n",
    "vx.plot(label='vx')\n",
    "vy.plot(label='vy')\n",
    "\n",
    "decorate(xlabel='Time (s)',\n",
    "         ylabel='Velocity (m/s)')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The x velocity slows down due to drag.\n",
    "\n",
    "The y velocity drops quickly while drag and gravity are in the same direction, then more slowly after the ball starts to fall.\n",
    "\n",
    "Another way to visualize the results is to plot y versus x.  The result is the trajectory of the ball through its plane of motion."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 40,
   "metadata": {},
   "outputs": [],
   "source": [
    "def plot_trajectory(results):\n",
    "    xs = results.R.extract('x')\n",
    "    ys = results.R.extract('y')\n",
    "    plot(xs, ys, color='C2', label='trajectory')\n",
    "\n",
    "    decorate(xlabel='x position (m)',\n",
    "             ylabel='y position (m)')\n",
    "\n",
    "plot_trajectory(results)\n",
    "savefig('figs/chap22-fig02.pdf')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Animation\n",
    "\n",
    "One of the best ways to visualize the results of a physical model is animation.  If there are problems with the model, animation can make them apparent.\n",
    "\n",
    "The ModSimPy library provides `animate`, which takes as parameters a `TimeSeries` and a draw function.\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The draw function should take as parameters a `State` object and the time.  It should draw a single frame of the animation.\n",
    "\n",
    "Inside the draw function, you almost always have to call `set_xlim` and `set_ylim`.  Otherwise `matplotlib` auto-scales the axes, which is usually not what you want."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 41,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = results.R.extract('x')\n",
    "ys = results.R.extract('y')\n",
    "\n",
    "def draw_func(state, t):\n",
    "    set_xlim(xs)\n",
    "    set_ylim(ys)\n",
    "    x, y = state.R\n",
    "    plot(x, y, 'bo')\n",
    "    decorate(xlabel='x position (m)',\n",
    "             ylabel='y position (m)')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 42,
   "metadata": {},
   "outputs": [],
   "source": [
    "animate(results, draw_func)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Delete the lines that set the x and y axes (or [comment them out](https://en.wiktionary.org/wiki/comment_out)) and see what the animation does."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Under the hood\n",
    "\n",
    "`Vector` is a function that returns a `ModSimVector` object."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 43,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = Vector(3, 4)\n",
    "type(V)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A `ModSimVector` is a specialized kind of Pint `Quantity`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 44,
   "metadata": {},
   "outputs": [],
   "source": [
    "isinstance(V, Quantity)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "There's one gotcha you might run into with Vectors and Quantities.  If you multiply a `ModSimVector` and a `Quantity`, you get a `ModSimVector`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 45,
   "metadata": {},
   "outputs": [],
   "source": [
    "V1 = V * m"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 46,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(V1)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "But if you multiply a `Quantity` and a `Vector`, you get a `Quantity`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 47,
   "metadata": {},
   "outputs": [],
   "source": [
    "V2 = m * V"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 48,
   "metadata": {},
   "outputs": [],
   "source": [
    "type(V2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a `ModSimVector` you can get the coordinates using dot notation, as well as `mag`, `mag2`, and `angle`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 49,
   "metadata": {},
   "outputs": [],
   "source": [
    "V1.x, V1.y, V1.mag, V1.angle"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "With a `Quantity`, you can't.  But you can use indexing to get the coordinates:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 50,
   "metadata": {},
   "outputs": [],
   "source": [
    "V2[0], V2[1]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And you can use vector functions to get the magnitude and angle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 51,
   "metadata": {},
   "outputs": [],
   "source": [
    "vector_mag(V2), vector_angle(V2)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "And often you can avoid the whole issue by doing the multiplication with the `ModSimVector` on the left."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** Run the simulation with and without air resistance.  How wrong would we be if we ignored drag?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 52,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hint\n",
    "\n",
    "system_no_drag = System(system, C_d=0)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 53,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 54,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 55,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 56,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 57,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** The baseball stadium in Denver, Colorado is 1,580 meters above sea level, where the density of air is about 1.0 kg / meter$^3$.  How much farther would a ball hit with the same velocity and launch angle travel?"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 58,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Hint\n",
    "\n",
    "system2 = System(system, rho=1.0*kg/m**3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 59,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 60,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**Exercise:** The model so far is based on the assumption that coefficient of drag does not depend on velocity, but in reality it does.  The following figure, from Adair, [*The Physics of Baseball*](https://books.google.com/books/about/The_Physics_of_Baseball.html?id=4xE4Ngpk_2EC), shows coefficient of drag as a function of velocity.\n",
    "\n",
    "<img src=\"data/baseball_drag.png\" width=\"400\">\n",
    "\n",
    "\n",
    "I used [an online graph digitizer](https://automeris.io/WebPlotDigitizer/) to extract the data and save it in a CSV file.  Here's how we can read it:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 61,
   "metadata": {},
   "outputs": [],
   "source": [
    "baseball_drag = pd.read_csv('data/baseball_drag.csv')\n",
    "mph = Quantity(baseball_drag['Velocity in mph'], UNITS.mph)\n",
    "mps = mph.to(m/s)\n",
    "baseball_drag.index = magnitude(mps)\n",
    "baseball_drag.index.name = 'Velocity in meters per second'\n",
    "baseball_drag"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Modify the model to include the dependence of `C_d` on velocity, and see how much it affects the results.  Hint: use `interpolate`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 62,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 63,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 64,
   "metadata": {},
   "outputs": [],
   "source": [
    "C_d = drag_interp(43 * m / s)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 65,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 66,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Solution goes here"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}
